Optimización heurística

Autores/as

Ronald Gabriel Palencia

Junior Antonio Muñoz Henao

Optimzación Heurìstica

Introdución

La optimización heurística engloba un conjunto de métodos destinados a encontrar soluciones aproximadas a problemas de optimización que son demasiado complejos para los enfoques analíticos tradicionales, como aquellos basados en el cálculo de gradientes. A diferencia de los métodos determinísticos, que procuran soluciones exactas bajo supuestos matemáticos estrictos sobre las funciones objetivo, los enfoques heurísticos exploran el espacio de soluciones de manera estratégica para encontrar soluciones satisfactorias sin garantizar necesariamente el óptimo global. Estos enfoques resultan particularmente valiosos en situaciones donde la función objetivo es discontinua, no diferenciable, altamente compleja, o cuando el espacio de búsqueda es tan amplio que hace inviable una exploración exhaustiva.

Dentro de los métodos heurísticos mas destacados y empleados se hallan los algoritmos evolutivos y las colonias de hormigas, ambos inspirados en fenómenos biológicos.

Algoritmos Evolutivos

Los algoritmos evolutivos son inspirados por el proceso de evolución biológica y selección natural, operando mediante un ciclo de selección, cruzamiento, y mutación en una población de soluciones candidatas (Mitchell, 1998). Cada solución es evaluada a través de una función de aptitud, determinando su probabilidad de ser seleccionada para la generación de nuevas soluciones. Este ciclo permite que, con el tiempo, las soluciones mas aptas predominen, dirigiendo al algoritmo hacia una solución óptima o satisfactoria.

Colonias de Hormigas

Por su parte, las colonias de hormigas simulan el comportamiento colectivo de las hormigas en su búsqueda de alimentos, donde cada hormiga explora el espacio de soluciones siguiendo caminos marcados por feromonas, sustancias químicas que permiten la comunicación entre ellas (Dorigo, Maniezzo, & Colorni, 1996). Las rutas con concentraciones mas altas de feromonas atraen a mas hormigas, favoreciendo la exploración de soluciones prometedoras a través de un proceso de retroalimentación positiva.

Aplicación de Métodos Heurísticos

La aplicación de métodos heurísticos para la optimización de funciones de prueba implica varios pasos clave, comenzando por la definición precisa de la función de prueba que representa el problema de optimización. La elección del método heurístico adecuado depende de las características del problema, incluyendo la naturaleza de la función objetivo y las peculiaridades del espacio de búsqueda. La configuración de los parámetros del algoritmo, como el tamaño de la población en los algoritmos evolutivos o el número de hormigas en las colonias de hormigas, es fundamental para la efectividad del método. Posteriormente, se ejecuta el algoritmo y se monitorea su progreso, realizando ajustes en los parámetros según sea necesario para optimizar los resultados.

El empleo de métodos heurísticos en la optimización presenta una alternativa robusta y adaptable para abordar problemas complejos, ofreciendo soluciones aproximadas cuando los enfoques convencionales resultan inadecuados o ineficaces. La inspiración en procesos naturales no solo proporciona una perspectiva innovadora para la resolución de problemas de optimización sino que también abre nuevas posibilidades de investigación y aplicación en áreas como ingeniería, economía, y ciencia de datos.

Parte 1: optimización numérica

Considere las siguientes funciones de prueba:

Función de Rosenbrock

La Función de Rosenbrock es conocida por su valle en forma de banana que contiene el mínimo global. Se define como:

\[ f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2 \]

con \(x_i \in [-2.048, 2.048]\), \(i = 1, 2\). Alcanza su valor mínimo en \(x_1 = 1\) y \(x_2 = 1\).

Función de Rastrigin

La Función de Rastrigin es una función no convexa utilizada como problema de prueba de rendimiento para algoritmos de optimización debido a su gran cantidad de mínimos locales. Se define como:

\[ f(x_1, x_2) = 20 + \sum_{i=1}^{2}(x_i^2 - 10 \cos(2\pi x_i)) \]

con \(x_i \in [-5.12, 5.12]\), \(i = 1, 2\). Alcanza su valor mínimo en \(x_1 = 0\) y \(x_2 = 0\).

Función de Schwefel

La Función de Schwefel es conocida por su complejidad y su paisaje de búsqueda desafiante. Se define como:

\[ f(x_1, x_2) = 418.9829 \times 2 - \sum_{i=1}^{2}x_i \sin(\sqrt{|x_i|}) \]

con \(x_i \in [-500, 500]\), \(i = 1, 2\). Alcanza su valor mínimo en \(x_1 = 420.9687\) y \(x_2 = 420.9687\).

Función de Griewank

La Función de Griewank a menudo se utiliza para probar la eficacia de algoritmos de optimización en tratar con oscilaciones grandes y numerosos óptimos locales. Se define como:

\[ f(x_1, x_2) = \frac{1}{4000}\sum_{i=1}^{2}x_i^2 - \prod_{i=1}^{2}\cos(\frac{x_i}{\sqrt{i}}) + 1 \]

con \(x_i \in [-600, 600]\), \(i = 1, 2\). Alcanza su valor mínimo en \(x_1 = 0\) y \(x_2 = 0\).

Función Goldstein-Price

La Función Goldstein-Price es conocida por tener un valle estrecho que conduce al mínimo global. Se define como:

\[ \begin{align*} f(x_1, x_2) = & [1 + (x_1 + x_2 + 1)^2(19 - 14x_1 + 3x_1^2 - 14x_2 + 6x_1x_2 + 3x_2^2)]\\ & \times [30 + (2x_1 - 3x_2)^2(18 - 32x_1 + 12x_1^2 + 48x_2 - 36x_1x_2 + 27x_2^2)] \end{align*} \]

con \(x_i \in [-2, 2]\), \(i = 1, 2\). Alcanza su valor mínimo en \(x_1 = 0\) y \(x_2 = -1\).

Función de las seis jorobas de camello (Six-Hump Camel Back)

La Función de las seis jorobas de camello es una función de prueba que presenta dos mínimos globales. Se define como:

\[ f(x_1, x_2) = (4 - 2.1x_1^2 + \frac{x_1^4}{3})x_1^2 + x_1x_2 + (-4 + 4x_2^2)x_2^2 \]

con (x_1 ) y (x_2 ). Alcanza su valor mínimo en (x_1 = -0.0898) y (x_2 = 0.7126), y también en (x_1 = 0.0898) y (x_2 = 0.7126).

Teniendo en cuanta las anteriores funciones realizar los siguientes items:

  1. Escoja dos funciones de prueba

  2. Optimice las funciones en dos y tres dimensiones usando un método de descenso por gradiente con condición inicial aleatoria

  3. Optimice las funciones en dos y tres dimensiones usando: algoritmos evolutivos, optimización de partículas y evolución diferencial

  4. Represente con un gif animado o un video el proceso de optimización de descenso por gradiente y el proceso usando el método heurístico.

Código 
#!pip install pygad
#!pip install numpy

import numpy as np

import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML;
import pygad
rc('animation', html='html5');

1. Funciones de Prueba

La funciones de prueba escogidas para probar los métodos de optimización serán las siguientes:

  • Función de Rastrigin

  • Función de Schwefel

Código 
# Se implementa la función de Rastrigin vectorizada para arreglos numpy
def Rastrigin(X):
  Y = 10 * len(X) + np.sum(X ** 2 - 10 * np.cos(2 * np.pi * X))
  return Y

# Se implementa la función de Schwefel vectorizada para arreglos numpy
def Schwefel(X):
    Y = np.sum(X * np.sin(np.sqrt(np.abs(X))))
    return(Y)

# Visualización de la función de Rastrigin en 2D
ncols = 100
nrows = 100
X1 = np.linspace(-5.12, 5.12, ncols)
Y1 = np.linspace(-5.12, 5.12, nrows)
X1, Y1 = np.meshgrid(X1, Y1)

Z1 = [Rastrigin(np.array([X1[i,j], Y1[i,j]])) for i in range(nrows) for j in range(ncols)]
Z1 = np.array(Z1).reshape([nrows,ncols])

fig, (ax1,ax2) = plt.subplots(1,2)
fig.set_size_inches(13,5)
contorno = ax1.contourf(X1,Y1,Z1)
ax1.set_title("Función de Rastrigin")
fig.colorbar(contorno, ax=ax1)
fig.subplots_adjust(wspace=1.5)

# Visualización de la función de Schwefel en 2D
X2 = np.linspace(-512, 512, ncols*10)
Y2 = np.linspace(-512, 512, nrows*10)
X2, Y2 = np.meshgrid(X2, Y2)
Z2 = [Schwefel(np.array([X2[i,j], Y2[i,j]])) for i in range(nrows*10) for j in range(ncols*10)]
Z2 = np.array(Z2).reshape([nrows*10,ncols*10])
contorno2 = ax2.contourf(X2,Y2,Z2)
fig.colorbar(contorno2, ax=ax2)
ax2.set_title("Función de Schewfel")


plt.tight_layout()
plt.show()

Función de Rastrigin (izquierda):

La gráfica muestra una superficie repleta de picos y valles regulares y simétricos, característicos de la Función de Rastrigin. Esta función es conocida por su gran número de mínimos locales que dificultan encontrar el mínimo global. Los colores mas oscuros representan valores mas bajos de la función, y los mas claros, valores mas altos. El patrón periódico de los mínimos locales es claramente visible y sugiere la complejidad de utilizar métodos de optimización que no se basan en gradientes para encontrar el óptimo global, que se sitúa en el centro de la gráfica donde se observa un mínimo profundo y distinto.

Función de Schwefel (derecha):

La gráfica para la Función de Schwefel muestra una estructura mas compleja con una serie de anillos concéntricos que representan los niveles de la función. La Función de Schwefel es conocida por tener un mínimo global pronunciado, pero rodeado de zonas que llevan a mínimos locales subóptimos, lo cual puede engañar a los algoritmos de búsqueda. Los valores mas bajos en esta función se indican con un color oscuro, ubicados cerca del centro de la gráfica, que es donde se encuentra el mínimo global. La amplia gama de colores refleja la amplia variedad de valores que la función puede tomar, lo que indica un paisaje de optimización desafiante con múltiples mínimos locales.

Visualización en 3D

Código 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Definición de las funciones
def Rastrigin(X):
    return 10 * len(X) + np.sum(X ** 2 - 10 * np.cos(2 * np.pi * X))

def Schwefel(X):
    return 418.9829 * len(X) - np.sum(X * np.sin(np.sqrt(np.abs(X))))

# Creación de la malla para graficar
ncols, nrows = 100, 100
X1 = np.linspace(-5.12, 5.12, ncols)
Y1 = np.linspace(-5.12, 5.12, nrows)
X1, Y1 = np.meshgrid(X1, Y1)
Z1 = np.array([Rastrigin(np.array([x, y])) for x, y in zip(np.ravel(X1), np.ravel(Y1))])
Z1 = Z1.reshape(X1.shape)

# Graficar Rastrigin en 3D
fig = plt.figure(figsize=(14, 7))

# Rastrigin
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(X1, Y1, Z1, cmap='viridis')
ax1.set_title("Función de Rastrigin")
ax1.set_xlabel('X1')
ax1.set_ylabel('Y1')
ax1.set_zlabel('Z1')

# Ajuste para la función de Schwefel
ncols, nrows = 1000, 1000
X2 = np.linspace(-500, 500, ncols)
Y2 = np.linspace(-500, 500, nrows)
X2, Y2 = np.meshgrid(X2, Y2)
Z2 = np.array([Schwefel(np.array([x, y])) for x, y in zip(np.ravel(X2), np.ravel(Y2))])
Z2 = Z2.reshape(X2.shape)

# Schwefel
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(X2, Y2, Z2, cmap='viridis')
ax2.set_title("Función de Schwefel")
ax2.set_xlabel('X2')
ax2.set_ylabel('Y2')
ax2.set_zlabel('Z2')

plt.tight_layout()
plt.show()

Las gráficas 3D muestran claramente la topología de las funciones de Rastrigin y Schwefel. Rastrigin tiene un paisaje ondulado con numerosos mínimos locales, mientras que Schwefel presenta un paisaje mas suave con un mínimo global prominente. Estas visualizaciones resaltan la complejidad inherente a cada función

2. Optimización mediante descenso por gradiente

A continuación, presentamos las implementaciones de las funciones necesarias para calcular el gradiente. Ademas, introducimos el método de descenso por gradiente, una técnica de optimización que utilizaremos para hallar los mínimos de las funciones objetivo. Este enfoque iterativo nos permitirá afinar progresivamente nuestras estimaciones hacia el punto de óptimo global.

Código 
# La siguiente función implementa el cálculo del gradiente usando 'list comprehension'
def num_grad(x, fun, h=0.01):
    return np.array([(fun(x+e*h) - fun(x-e*h)) / (2*h) for e in np.eye(len(x))])

# La siguiente función implmenta la optimización mediante descenso por gradiente
def optimizador_mult_numdev(f, x0, eta=0.01, tol=1e-6, max_iter=1e3):
    x = np.array(x0)
    iter = 0
    while np.linalg.norm(num_grad(x, f, h=0.01)) > tol and iter < max_iter:
        x = x - eta*num_grad(x, f, h=0.01)
        iter += 1
    return x

2.1 Función de Rastrigin

Generación del gráfico

Se procede a definir las características del gráfico que representará todo.

::: {#059e9496 .cell execution_count=5} ``` {.python .cell-code} x0 = np.array([4np.random.random()-1,4np.random.random()-1]) xs = [x0]

for i in range(25): x_new = optimizador_mult_numdev(Rastrigin, x0) xs.append(x_new) x0 = x_new ``` :::

Código 
# Particionamiento del rango de cada variable
ncols = 100
nrows = 100
X = np.linspace(-3.12, 3.12, ncols)  #Posibles valores de X
Y = np.linspace(-3.12, 3.12, nrows)  #Posibles valores de Y
X, Y = np.meshgrid(X, Y)  #Definición del plano con posibles valores

# Evaluación de la función de Rastrigin
Z = [Rastrigin(np.array([X[i,j], Y[i,j]])) for i in range(nrows) for j in range(ncols)]
Z = np.array(Z).reshape([nrows,ncols])

fig, ax = plt.subplots()
fig.set_size_inches(10,7)
contorno = ax.contour(X,Y,Z)  #Curvas de nivel de la función de Rastrigin
contornof = ax.contourf(X,Y,Z) #Rellenando las curvas de nivel
line, = ax.plot(xs[-1][0], xs[-1][1], 'xr--',mec='b', markersize=10) #Creando la línea que contendrá los puntos
              # de la animación de Optimización
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Función de Rastrigin')
plt.colorbar(contornof) # Se muestra la barra lateral con la escala de valores
          # ¡¡Éste método depende de la creación de un objeto "contourf" !!
plt.show()

En el gráfico anterior, se destaca el punto que representa el mínimo encontrado a través del método de descenso por gradiente. A continuación, presentaremos un GIF animado que muestra el conjunto de soluciones halladas en cada iteración del proceso hasta alcanzar la solución óptima. Esta visualización animada brinda una perspectiva clara y dinámica de cómo el algoritmo de descenso por gradiente avanza iterativamente a través del espacio de soluciones, acercándose paso a paso al punto de mínimo global de la función de Rastrigin. Este enfoque no solo facilita la comprensión del proceso de optimización, sino que también ilustra de manera efectiva la eficacia del algoritmo en la búsqueda de soluciones óptimas.

Código 
from matplotlib import rc
rc('animation', html='jshtml')

# Definición de la función graficadora
def graficar(frame):
    x, y = zip(*xs[:frame+1])
    line.set_data(x, y)
    return line,


# Animación final
ani = animation.FuncAnimation(fig, graficar, frames=len(xs), interval=300)
ani

Optimización en 3D

La siguiente gráfica ilustra el punto mínimo de la función de Rastrigin en tres dimensiones. Es notable que este mínimo está destacado en color rojo.

Código 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def Rastrigin(x):
    A = 10
    return A * len(x) + sum([(xi ** 2 - A * np.cos(2 * np.pi * xi)) for xi in x])

def num_grad(x, fun, h=0.001):
    return np.array([(fun(x + e * h) - fun(x - e * h)) / (2 * h) for e in np.eye(len(x))])

def optimizador_mult_numdev(f, x0, eta=0.0001, tol=1e-6, max_iter=1e4):
    x = np.array(x0)
    iter = 0
    while np.linalg.norm(num_grad(x, f)) > tol and iter < max_iter:
        x = x - eta * num_grad(x, f)
        iter += 1
    return x

ncols, nrows = 100, 100
X1 = np.linspace(-5.12, 5.12, ncols)
Y1 = np.linspace(-5.12, 5.12, nrows)
X1, Y1 = np.meshgrid(X1, Y1)
Z1 = np.array([Rastrigin(np.array([x, y])) for x, y in zip(np.ravel(X1), np.ravel(Y1))])
Z1 = Z1.reshape(X1.shape)

# Ajustar el punto inicial si es necesario
xmin = optimizador_mult_numdev(Rastrigin, [-3.0, 3.0], eta=0.01, max_iter=10000)  
zmin = Rastrigin(xmin)

# Graficar
fig = plt.figure(figsize=(14, 7))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X1, Y1, Z1, cmap='viridis', edgecolor='none', alpha=0.7)
ax.scatter(xmin[0], xmin[1], zmin, color='r', s=50, label=f'Punto mínimo encontrado en {xmin} con valor {zmin}')
ax.set_title("Función de Rastrigin con punto mínimo")
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()

plt.show()

El comportamiento del optimizador de descenso por gradiente, partiendo de un punto inicial \((-3, 3)\) y moviéndose hacia aproximadamente \((-2.96, 2.96)\) con un valor de función de 41, evidencia las limitaciones al enfrentar la compleja función de Rastrigin. Esta mejora, sin alcanzar el mínimo global, subraya las dificultades debido a la naturaleza multimodal de Rastrigin, caracterizada por abundantes mínimos locales. El resultado demuestra la capacidad del optimizador para progresar en un espacio desafiante, utilizando información local del gradiente en su búsqueda de optimización.

2.2 Función de Schwefel

texto del vínculoCon las funciones de derivadas parciales y

de optimización ya implementadas, sólo resta recalcular los datos para la función de Schwefel y generar el gráfico.

Código 
def optimizador_mult_numdev(f, x0, eta=0.0001, tol=1e-6, max_iter=1e4):
    x = np.array(x0)
    iter = 0
    while np.linalg.norm(num_grad(x, f)) > tol and iter < max_iter:
        x = x - eta * num_grad(x, f, h = 1e-4)
        iter += 1
    return x
np.random.seed(42)
# Implementación del método de descenso por gradiente
x0 = np.array([300*(2*np.random.random()-1),300*(2*np.random.random()-1)]) # generación de un punto aleatorio
#print(x0)
#xs = [x0]
for i in range(25):
    x_new = optimizador_mult_numdev(Schwefel, x0)
    xs.append(x_new)
    x0 = x_new

print(xs)
[array([1.93429234, 1.02608422]), array([1.78610916, 1.18287674]), array([1.81950437, 0.66108662]), array([1.81997582, 0.86847643]), array([1.82014478, 1.22197859]), array([1.80657218, 0.92227481]), array([1.75906982, 0.87537773]), array([1.79859281, 1.32599757]), array([1.8329174 , 1.38230191]), array([1.81291812, 0.89576452]), array([1.83710258, 0.69712025]), array([1.83548921, 1.27741073]), array([1.83462804, 0.94766826]), array([1.81768718, 1.28689829]), array([1.76213092, 0.87265638]), array([1.82130468, 1.36516073]), array([1.83480425, 0.58115991]), array([1.83301351, 0.57808453]), array([1.82518858, 1.32517258]), array([1.83389648, 0.65215772]), array([1.78572404, 0.72616781]), array([1.79105771, 0.87163929]), array([1.8373511 , 1.28240507]), array([1.80337715, 0.637759  ]), array([1.8083492 , 0.96791667]), array([1.81115952, 1.25181321]), array([-78.07857198, 263.25135889]), array([-81.52866165, 255.45937899]), array([-85.62282263, 247.49152593]), array([-90.24330614, 239.84467053]), array([-95.14570071, 232.93049206]), array([-100.00989018,  226.9825909 ]), array([-104.5362827 ,  222.05457484]), array([-108.52476498,  218.07659714]), array([-111.89468087,  214.91950269]), array([-114.65771164,  212.44010227]), array([-116.87736877,  210.50510049]), array([-118.63667915,  209.00041641]), array([-120.01899781,  207.8327029 ]), array([-121.09900602,  206.92746076]), array([-121.93974797,  206.22606055]), array([-122.59267692,  205.68272204]), array([-123.09895385,  205.26185159]), array([-123.49110628,  204.93583873]), array([-123.79464381,  204.68329011]), array([-124.0294754 ,  204.48763811]), array([-124.21109015,  204.33605443]), array([-124.35151351,  204.21860614]), array([-124.46006878,  204.1276015 ]), array([-124.54397738,  204.05708364]), array([-124.608829  ,  204.00243867])]
Código 
# Particionamiento del rango de cada variable
ncols = 1000
nrows = 1000
X = np.linspace(-512, 512, ncols)  #Posibles valores de X
Y = np.linspace(-512, 512, nrows)  #Posibles valores de Y
X, Y = np.meshgrid(X, Y)  #Definición del plano con posibles valores

# Evaluación de la función de Rastrigin
Z = [Schwefel(np.array([X[i,j], Y[i,j]])) for i in range(nrows) for j in range(ncols)]
Z = np.array(Z).reshape([nrows,ncols])

fig, ax = plt.subplots()
fig.set_size_inches(10,7)
contorno = ax.contour(X,Y,Z)  #Curvas de nivel de la función de Rastrigin
contornof = ax.contourf(X,Y,Z) #Rellenando las curvas de nivel
line, = ax.plot(xs[-1][0], xs[1][1], 'xr--',mec='b') #Creando la línea que contendrá los puntos
              # de la animación de Optimización
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Función de Schwefel')
plt.colorbar(contornof) # Se muestra la barra lateral con la escala de valores
          # ¡¡Éste método depende de la creación de un objeto "contourf" !!
plt.show()

Durante la ejecución del optimizador de descenso por gradiente aplicado a la función de Schwefel, el punto de inicio aleatorio seleccionado fue \((-75.28, 270.43)\). El optimizador convergió hacia el punto \((-65.56, 302.46)\), que está significativamente alejado del mínimo global conocido de la función. Esta dificultad para alcanzar el mínimo global puede atribuirse a la presencia de numerosos mínimos locales en la función de Schwefel, un desafío común para los métodos de optimización basados en el gradiente. Ademas, se realizó una exploración variando el valor de \(\eta\), observándose que la reducción de este parámetro resultaba en una ejecución mas lenta del algoritmo. En una sección posterior, se ilustrará este proceso de optimización mediante un GIF, que ofrecerá una visualización de la trayectoria seguida por el optimizador en su búsqueda de los valores mínimos encontrados.

Código 
# Implementación del método de descenso por gradiente
x0 = np.array([300*(2*np.random.random()-1),300*(2*np.random.random()-1)]) # generación de un punto aleatorio
xs = [x0]
for i in range(25):
    x_new = optimizador_mult_numdev(Schwefel, x0)
    xs.append(x_new)
    x0 = x_new

# Definición de la función graficadora
def graficar(frame):
    x, y = zip(*xs[:frame+1])
    line.set_data(x, y)
    return line,

# Animación final
ani = animation.FuncAnimation(fig, graficar, frames=len(xs), interval=300)
ani

Optimización en 3D

La optimización en tres dimensiones de la función de Schwefel no consigue ubicar el mínimo global, sin embargo, evidencia un avance significativo hacia este. La gráfica subsiguiente ilustra la secuencia de puntos descubiertos durante el proceso de optimización, lo cual demuestra el progresivo acercamiento hacia la región de mínimo valor en la función objetivo.

Código 
#!pip install plotly
Código 
np.random.seed(42)
x0 = np.array([300*(2*np.random.random()-1) for _ in range(3)])  # Punto inicial en 3D
print(f"Punto inicial: {x0}")

xs = [x0]  # Almacenará todos los puntos por los que pasa el optimizador
for i in range(25):
    x_new = optimizador_mult_numdev(Schwefel, x0, eta=0.0001, tol=1e-6, max_iter=10000)
    xs.append(x_new)
    x0 = x_new

xs = np.array(xs)  # Convertir a un array de Numpy para facilitar la manipulación
Punto inicial: [-75.27592869 270.42858385 139.19636509]
Código 
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# Trayectoria del optimizador
ax.plot(xs[:,0], xs[:,1], xs[:,2], marker='o', linestyle='-', color='r', label='Trayectoria del Optimizador')

ax.set_title('Trayectoria del Optimizador de Descenso por Gradiente en la Función de Schwefel')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()

plt.show()

Mediante la visualización interactiva en tres dimensiones, se confirma que, al igual que en el caso bidimensional, el optimizador no alcanza el mínimo global de la función de Schwefel. Este resultado es consecuencia de la abundancia de mínimos locales que caracterizan a esta función. No obstante, es interesante observar que, pese a la selección aleatoria del punto inicial, el optimizador no queda atrapado en ninguno de estos mínimos locales, lo cual sugiere una dinámica de búsqueda que consigue evitar estancamientos prematuros en la exploración del espacio de soluciones.

Código 
import plotly.graph_objs as go
import numpy as np

# Definir la función de Schwefel para tres dimensiones
def Schwefel(x):
    return 418.9829 * len(x) - sum([xi * np.sin(np.sqrt(abs(xi))) for xi in x])

# Generar datos para la superficie de la función de Schwefel
x = np.linspace(-500, 500, 100)
y = np.linspace(-500, 500, 100)
X, Y = np.meshgrid(x, y)
Z = np.array([Schwefel([X[i][j], Y[i][j], 0]) for i in range(X.shape[0]) for j in range(X.shape[1])])
Z = Z.reshape(X.shape)

# Generar el punto optimizado (aquí necesitas tus valores optimizados)
x_opt = [420.9687, 420.9687, 0]  # Reemplaza con el resultado de tu optimización
z_opt = Schwefel(x_opt)

# Crear la figura
fig = go.Figure(data=[
    go.Surface(z=Z, x=X, y=Y, colorscale='Viridis'),
    go.Scatter3d(x = xs[:,0], y = xs[:,1], z = xs[:,2], mode='markers', marker=dict(size=5, color='red'))
])

# Actualizar el layout de la figura
fig.update_layout(title='Función de Schwefel con Punto Optimizado en 3D', autosize=False,
                  width=700, height=700,
                  margin=dict(l=65, r=50, b=65, t=90))

# Mostrar la figura
fig.show()

3. Optimización mediante algoritmos evolutivos

3.1 Función de Rastrigin

- Optimización de la función de Rastrigin

Se configura un algoritmo genético mediante la biblioteca PyGAD para la optimización de la función de Rastrigin. La configuración se detalla de la siguiente manera:

  • El algoritmo iterará a través de 60 generaciones.

  • En cada generación se cruzarán dos individuos.

  • La función objetivo que guiará la adaptación de los individuos es la propia función de Rastrigin.

  • Cada generación estará compuesta por una población de 9 individuos.

  • La selección para la reproducción se basará en la elección de los individuos mas aptos.

  • La mutación de los individuos seguirá una estrategia aleatoria.

  • Esta configuración establece los parámetros operativos fundamentales del algoritmo genético y se orienta a la optimización efectiva de la compleja función de Rastrigin.

Código 
#!pip install pygad==3.3.0

Una vez creada la instancia de nuestro algoritmo, ejecutamos la optimización mediante el algoritmo evolutivo y mostramos la evolución del ajuste generación tras generación

ga_instance_rastrigin = pygad.GA(num_generations=1, # We’ll manually loop for more generations num_parents_mating=2, fitness_func=Rastrigin_ga, sol_per_pop=10, num_genes=3, init_range_low=-5.12, init_range_high=5.12)

Código 
import pygad

# Define your fitness function
def Rastrigin_ga(solution, solution_idx):
    X = solution
    Y = 10 * len(X) + sum([(x ** 2 - 10 * np.cos(2 * np.pi * x)) for x in X])
    fitness = -Y
    return fitness

# Initialize the GA
ga_instance_rastrigin = pygad.GA(num_generations=1,  # We'll manually loop for more generations
                                 num_parents_mating=2,
                                 fitness_func=Rastrigin_ga,
                                 sol_per_pop=10,
                                 num_genes=3,
                                 init_range_low=-5.12,
                                 init_range_high=5.12)

fitness_history = []

# Manually run the GA for multiple generations
for generation in range(60):
    ga_instance_rastrigin.run()
    best_solution, best_solution_fitness, best_solution_idx = ga_instance_rastrigin.best_solution()
    fitness_history.append(best_solution_fitness)
    # print("Generation : ", generation + 1)
    # print("Fitness of the best solution :", best_solution_fitness)

# You can plot the fitness history using Matplotlib if desired.
import matplotlib.pyplot as plt
plt.plot(fitness_history)
plt.title('Fitness Evolution Over Generations')
plt.xlabel('Generation')
plt.ylabel('Fitness')
plt.show()

En la gráfica anterior, se puede observar la trayectoria de optimización del algoritmo genético PyGAD en función del número de generaciones. Se aprecia una marcada mejora en la aptitud durante las etapas iniciales, con un rápido incremento hacia valores mas óptimos. Este progreso temprano sugiere que el algoritmo es eficiente en la identificación de soluciones prometedoras desde el principio. Tras esta fase de mejoras significativas, la curva se nivela, lo cual indica que el algoritmo ha llegado a un punto de convergencia y se están haciendo menos descubrimientos de nuevas soluciones óptimas a medida que avanza el número de generaciones. La estabilidad de la aptitud hacia el final de las iteraciones sugiere que se ha alcanzado un punto cercano al óptimo, donde las subsiguientes evoluciones ofrecen ganancias marginales en el desempeño o aptitud del modelo.

- Muestra de gráficos

La gráfica siguiente exhibe el punto óptimo encontrado por el algoritmo genético, marcado con una ‘X’, situado en una zona de bajo valor, que indica un rendimiento superior en la tarea de optimización en comparación con el método de descenso por gradiente. Este resultado resalta la habilidad del algoritmo genético para evitar los mínimos locales y dirigirse hacia una solución mas cercana al mínimo global, lo que demuestra su eficacia en la exploración del espacio de soluciones y su ventaja sobre métodos basados en el gradiente en escenarios de optimización complejos.

Código 
import pygad
import numpy as np

def Rastrigin_ga(solution, solution_idx):
    X = solution
    Y = 10 * len(X) + np.sum(X ** 2 - 10 * np.cos(2 * np.pi * X))
    return -Y  # We will maximize this, hence the negative sign

# Define the GA
ga_instance = pygad.GA(num_generations=60,
                       num_parents_mating=2,
                       fitness_func=Rastrigin_ga,
                       sol_per_pop=10,
                       num_genes=2,
                       init_range_low=-5.12,
                       init_range_high=5.12)

# Run the GA
ga_instance.run()

# Collect the best solution after the final generation
best_solution = ga_instance.best_solution()[0]

print(best_solution)
[-0.00578913 -0.00159222]
Código 
import matplotlib.pyplot as plt

# Prepare the grid and evaluate the function
ncols, nrows = 100, 100
X = np.linspace(-5.12, 5.12, ncols)
Y = np.linspace(-5.12, 5.12, nrows)
X, Y = np.meshgrid(X, Y)
Z = np.array([Rastrigin_ga(np.array([x, y]), None) for x, y in zip(X.ravel(), Y.ravel())]).reshape(X.shape)

# Plotting
fig, ax = plt.subplots()
fig.set_size_inches(10, 7)
contour = ax.contour(X, Y, Z, levels=50)
contourf = ax.contourf(X, Y, Z, levels=50)
plt.colorbar(contourf, ax=ax)

# Plot the best solution
ax.plot(best_solution[0], best_solution[1], 'rx', markersize=10, label='Best Solution')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Rastrigin Function with GA Optimization')
ax.legend()

plt.show()

3.2 Función de Schwefel

- Optimización de la función de Schwefel

Aquí se crea una instancia de nuestro algoritmo genético usando pygad para optimizar la función de Schwefel.

Especificaciones de los parámetros:

  • Se usa un número de generaciones de 60
  • Se emparejan 2 individuos en cada generación
  • La función de ajuste será nuestra función de prueba
  • Se tomarán 9 individuos en cada generación
  • El criterio de selección de reproducción : mejores individuos
  • Tipo de mutación: Aleatoria
Código 
#!pip install --upgrade pygad
Código 
import pygad

# Define your fitness function
def Schwefel_ga(X,solutionidx):
    Y = np.sum(X * np.sin(np.sqrt(np.abs(X))))
    fitness = -Y
    return(fitness)

# Initialize the GA
ga_instance_rastrigin = pygad.GA(num_generations=1,  # We'll manually loop for more generations
                                 num_parents_mating=2,
                                 fitness_func=Schwefel_ga,
                                 sol_per_pop=10,
                                 num_genes=3,
                                 init_range_low=-5.12,
                                 init_range_high=5.12)

fitness_history = []

# Manually run the GA for multiple generations
for generation in range(60):
    ga_instance_rastrigin.run()
    best_solution, best_solution_fitness, best_solution_idx = ga_instance_rastrigin.best_solution()
    fitness_history.append(best_solution_fitness)
    # print("Generation : ", generation + 1)
    # print("Fitness of the best solution :", best_solution_fitness)

# You can plot the fitness history using Matplotlib if desired.
import matplotlib.pyplot as plt
plt.plot(fitness_history)
plt.title('Fitness Evolution Over Generations')
plt.xlabel('Generation')
plt.ylabel('Fitness')
plt.show()

Tras configurar nuestra instancia del algoritmo genético, orientada específicamente a la optimización de la función de Schwefel, procedemos con la ejecución del proceso evolutivo. Durante esta fase, el algoritmo emprende una busqueda sistemática a través de las generaciones para optimizar la adaptación de las soluciones. A medida que avanza, capturamos y visualizamos el progreso de la aptitud de las soluciones, generación tras generación, ilustrando así la eficacia del algoritmo en su camino hacia la localización de un óptimo cercano al mínimo global. Este proceso nos permite observar de forma clara cómo el algoritmo genético mejora y refina las soluciones con cada iteración.

- Muestra de gráficos

Tras regenerar el gráfico y marcar el valor óptimo encontrado, se observa que este dista del verdadero mínimo global. Esta discrepancia subraya la complejidad inherente a la función de Schwefel, que, debido a su paisaje lleno de mínimos locales, plantea un desafío significativo para el algoritmo. No obstante, ajustando y experimentando con distintos parámetros dentro de la configuración de pygad.GA, es plausible que se mejore la capacidad del algoritmo para acercarse mas al mínimo global. Tales ajustes pueden incluir la modificación del tamaño de la población, la tasa de mutación, o la estrategia de selección, entre otros, lo cual podría optimizar la búsqueda y posiblemente conducir a resultados mas precisos.

Código 
import pygad
import numpy as np
def Schwefel_ga(X,solutionidx):
    Y = np.sum(X * np.sin(np.sqrt(np.abs(X))))
    fitness = -Y
    return(fitness)
# Define the GA
ga_instance = pygad.GA(num_generations=60,
                       num_parents_mating=2,
                       fitness_func=Schwefel_ga,
                       sol_per_pop=10,
                       num_genes=2,
                       init_range_low=-5.12,
                       init_range_high=5.12)

# Run the GA
ga_instance.run()

# Collect the best solution after the final generation
best_solution = ga_instance.best_solution()[0]
print(best_solution)
[-5.23584773 -5.24120212]
Código 
import matplotlib.pyplot as plt
import numpy as np

# Definir la función de fitness modificada para graficar
def Schwefel_ga(ga_instance, X, solution_idx):
    Y = np.sum(X * np.sin(np.sqrt(np.abs(X))))
    fitness = -Y
    return fitness

# Preparar la cuadrícula y evaluar la función
ncols, nrows = 100, 100
X = np.linspace(-5.12, 5.12, ncols)
Y = np.linspace(-5.12, 5.12, nrows)
X, Y = np.meshgrid(X, Y)

# Evaluar la función Schwefel en la cuadrícula, pasando `None` para `ga_instance` y un valor arbitrario para `solution_idx`
Z = np.array([Schwefel_ga(None, np.array([x, y]), 0) for x, y in zip(X.ravel(), Y.ravel())]).reshape(X.shape)

# Graficar
fig, ax = plt.subplots()
fig.set_size_inches(10, 7)
contour = ax.contour(X, Y, Z, levels=50)
contourf = ax.contourf(X, Y, Z, levels=50)
plt.colorbar(contourf, ax=ax)

# Graficar la mejor solución
ax.plot(best_solution[0], best_solution[1], 'rx', markersize=10, label='Best Solution')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Schwefel Function with GA Optimization')
ax.legend()

plt.show()